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Abstract 

We implement and investigate the numerical properties of a new family 
of integrators connecting both variants of the symplectic Euler schemes, 
and including an alternative to the classical symplectic mid-point scheme, 
with some additional terms. This family is derived from a new algorithm, 
introduced in a previous study, for generating symplectic integrators based 
on the concept of special symplectic manifold. The use of symplectic 
rotations and a particular type of projection keeps the whole procedure 
within the symplectic framework. 

We show that it is possible to define a set of parameters that control 
the additional terms providing a way of “tuning” these new symplectic 
schemes. We test the “tuned” symplectic integrators with the perturbed 
pendulum and we compare its behavior with an explicit SABA2 scheme 
for perturbed systems. Remarkably, for the given examples, the error 
in the energy integral can be reduced considerably. There is a natural 
geometrical explanation, sketched at the end of this paper. This is the 
subject of a parallel article where a finer analysis is performed. Numerical 
results obtained in this paper open a new point of view on symplectic 
integrators and Hamiltonian error. 

*The first author is supported by a grant from the Fondation du College de France under 
the research convention PU14150472. 
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1 Introduction 


A symplectic integrator for a Hamiltonian system is a numerical method which 
preserves the structure of the Hamiltonian vector field. Poincare discovered that 
the flow of a Hamiltonian system forms a one-parameter subgroup of canonical 
transformations. In modern language, we say that the Hamiltonian flow is 
symplectic. The standard procedure for simulating Hamiltonian dynamics is by 
discretizing the Hamiltonian flow, which consists in discretizing the evolution 
time and looking for symplectic transformations which map the state of the 
system between two adjacent elements of the discretized time, e.g. from time 
t n to time t n + 1 - 

It is well-known that one method for creating symplectic maps is based 
on generating functions, and it was already used by Poincare when looking 
for periodic orbits of second genus [30]. In fact, generating functions were an 
important ingredient of the invariant integral theory pioneered by Poincare and 
generalized by Cartan [4]. Since then, interest in generating functions remains 
very active from both the theoretical and the numerical point of view. Indeed 
for every symplectic transformation ip there corresponds (at least locally) a class 
of functions generating a Lagrangian submanifold, i.e., the graph of a 1-form 
symplectomorphic to the Liouville form. In addition, Hamilton-Jacobi theory 
connects this Lagrangian submanifold with another submanifold invariant under 
the transformation ip. 

Lagrangian submanifolds used to obtain suitable maps for symplectic inte¬ 
grators must have a very particular form. They must contain all the information 
concerning the source and the target (symplectic) variables. This information, 
encoded in the Lagrangian submanifold, contains the well-known fact that gen¬ 
erating functions for symplectic integrators must generate maps close to the 
identity. This criterion is not enough for determining whether the generating 
function associated to some symplectic map is suitable for obtaining a symplec¬ 
tic integrator. What is generally missing in the literature, is a geometrical ap¬ 
proach distilling the theory behind the different techniques within a unified point 
of view on symplectic integrators. The present paper, together with [18, 15] are 
contributions in this direction. Before going over the theory, we briefly review 
previous results on symplectic integrators and generating functions related to 
the present paper. 

The first article dealing with symplectic algorithms is attributed to De Voge- 
laere [38] in 1956. In 1983 Ruth [31] and Channell [5] made some progress with 
different techniques, in particular, Ruth pioneered explicit symplectic integra¬ 
tors and composition-splitting methods. Additional contributions were made 
independently by Menyuk [26] and Kang [19] in 1984. The same year, Kang 
Feng and his collaborators started a systematic study of symplectic integrators 
using generating functions [20, 13, 9]. His point of view was mostly algebraic, 
and based on Siegel’s article [33], reprinted some years later in book format [34]. 
Some geometrization was achieved by Ge and Marsden [13], Ge [10, 11], Ge and 
Dau-liu [12], although their numerical algorithms were based on Feng’s proce- 
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dure 1 . In 1990 Channell and Scovel [29] introduced symbolic computations to 
derive the integration formulas which arose in the procedure. Other important 
contributions were made by Sanz-Serna [32] who worked on symplecticity con¬ 
ditions for Runge-Kutta methods, and Miesbach and Pesch [27] who introduced 
some methods using Runge-Kutta techniques with generating functions. Many 
other authors have produced algorithms using generating functions, but the ge¬ 
ometrical construction remains the same, and their contributions diverge from 
our discussion. 

In a recent work [18], another link between generating functions for symplec- 
tic integrators and symplectic geometry has been studied, based on the concept 
of special symplectic manifold introduced by Tulczjyew in [40, 36]. Generating 
forms and functions in this framework were studied by Sniatycki and Tulczjyew 
[35] and Benenti [2, 3]. However, contributions from many other authors play 
a central role for the development and understanding of generating functions 
and their relationship with the Hamilton-Jacobi theory, such as Viterbo [37], 
Chaperon [6] Maslov [25], Hormander [14], Weinstein [39] among many others. 

In [18], the first author gives a strong argument for the construction of sym¬ 
plectic integrators based on the fact that solutions of the Hamilton-Jacobi equa¬ 
tion for integrable and autonomous systems belong to a Lagrangian submanifold 
of the phase space [1]. Starting from the classical approach, the product man¬ 
ifold of two copies of the phase space is created in [18]. Then a generalized 
generating function with 4 n variables is defined on the open ball of the product 
manifold centered at the initial condition (the source point of the map) and 
such that it contains the target point. The generating function is directly as¬ 
sociated to a primitive 1-form on the product manifold considered as a special 
symplectic manifold on the configuration space. The problem is translated from 
looking for the generating function, to looking for the 1-form, also known as 
the Liouvillian form [28]. Two different Lagrangian submanifolds arise, one de¬ 
fined by the generating function, solution of the Hamilton-Jacobi equation, and 
the other invariant under the flow of the Hamiltonian vector field. Applying 
an analogous argument to the Hamilton’s method of characteristics, the flow is 
searched in a transversal direction to the former submanifold. The suplemen- 
tary space, transversal to the tangent space of the first Lagrangian submanifold, 
is projected by the induced projection 2 , onto a 2 n subspace, where the original 
Hamiltonian system is finally evaluated. The projection induces a family of 
one step implicit symplectic integrators of generic order l 3 , closely related to 
those already studied by Kang and co-workers and recently revisited by Xue 
and Zanna [41]. 

The main difference is that our numerical schemes contain some additional 
terms which vanish when the stepsize goes to zero. In particular, the symmet- 

1 Most of Feng’s articles and some from his collaborators were recently edited in book 
format by M. Qin [21], 

2 The induced projection is the one we used to define the symplectic form on the product 
manifold by its pull-back. 

3 For some particular values of the parameters, we observe an increment in the order of 
convergence (see Figure 8). 


3 




ric integrators of the new family depend on 2 n parameters and contain, as a 
special case, the midpoint rule. The goal of the present article is to perform 
a numerical study investigating those additional terms and their influence on 
the accuracy and performance. It is a numerically oriented paper; for the geo¬ 
metrical arguments and theoretical development we refer the reader to [18, 16]. 
More theoretical results related to the numerical error and the claim by Ge and 
Marsen [13, 11] about the impossibility of constructing exact numerical sym- 
plectic mappings is addressed in [15]. In fact, Ge’s result [11] implies that our 
scheme should be exact up to numerical computer error and the residual error 
must be associated with the solution of the implicit equations. 

2 Hamiltonian systems and Symplectic integra¬ 
tors 

In what follows, we assume the reader is familiar with the terminology of differ¬ 
ential geometry and vector bundles. For an introduction, the reader is referred 
to [1, 23, 24], 

The approach presented in this work is based on the fact that Hamilto¬ 
nian mechanics relies on the geometrical properties present in the evolution of 
a mechanical system which accepts a Hamiltonian description. For this reason, 
in this section we use the standard notation in modern symplectic geometry, 
topology and fiber bundles. The goal is to give a formal geometrical framework 
to study symplectic integrators as isometries of a generic symplectic form on 
generic symplectic manifolds. In this way, the construction of our symplectic 
methods takes a more abstract and general point of view solving, or correcting, 
some misunderstanding arising when the analysis is restricted to symplectic vec¬ 
tor spaces, or cotangent bundles of linear spaces. In general, the starting point 
in symplectic integrator’s analysis is the identification of a cotangent bundle 
with a symplectic vector space by the isomorphisms 

T*(r) =° (R”)* © R" R 2n , 

where T*{ R ra ) is the cotangent bundle and (R n )* is the dual space of R”. How¬ 
ever, this point of view hides the geometrical background of generating functions 
for constructing symplectic maps. Let us start with the main definitions and 
results. 

A symplectic manifold is a 2n-dimensional manifold M equipped with a 
non-degenerated, skew-symmetric, closed 2-form u, such that at every point 
m € M, the tangent space to M at m, denoted T m M, has the structure of a 
symplectic vector space. One of the basic properties in symplectic geometry is 
given by Darboux’s theorem which states that any symplectic manifold is locally 
symplectomorphic to a symplectic vector space (V, cuo) with the canonical sym¬ 
plectic form wo(-,-) = (•, Jq-), where (•,•) is the canonical Euclidean structure 


4 




on V and Jq is the almost complex structure represented by the matrix: 


J 0 = 


0 „ 

In 



0 n , In £ M nxn (R). 


Jo is also known as the canonical symplectic matrix on V. Consequently, the 
tangent space to M at m, with its symplectic form w| m = w m and a suitable 
change of coordinates, can be completely described by the symplectic vector 
space with canonical symplectic coordinates (x, y) €V,yi = JoXi. Darboux’s 
theorem means that we systematically identify 


(T m M, Um ) = ( V ., w 0 ), Vm £ M, 


selecting the right change of symplectic coordinates to describe the dynamics 
on M by the canonical symplectic coordinates (x, y) e V. 

Remark 1 Unfortunately, Darboux’s theorem hides a very rich environment 
suitable for investigating symplectic maps by use of generating functions and 
Liouvillian forms. In practice, it locally identifies all the symplectic manifolds 
of the same dimension with the cotangent bundle. The information about the 
symplectomorphism which maps the symplectic manifold of interest to the cotan¬ 
gent bundle is lost when we apply Darboux’s theorem. To avoid this loss of in¬ 
formation we stay in the generic geometric framework of symplectic geometry 
by using special symplectic manifolds. 

In this geometrical framework, a Hamiltonian system (M, w, Xh) is a vector 
field A' = XH on the symplectic manifold (M, ui) such that 4 

ix = -dH, (1) 


for a differentiable function H : M —> R. 

A natural diffeomorphism b : TM —»• T*M between the tangent and the 
cotangent bundles of M is given by the contraction of the symplectic form with 
the vector field in the following way 


X ^ u(X, •)• 

The inverse of b is denoted by ft : T*M —> TM. Using ft, equation (1) is written 
in vector field form as 


A H = JVH, (2) 

where V is the standard gradient associated to the Euclidean structure. Ex¬ 
pression (2) is best suited for applications. The equations of evolution can be 
written as 


z = JV Z H( z), z € M. (3) 

4 Some authors write ix = dH instead of (1), but it depends on the definition of lj as 
2-form and the choice of the complex structure J. 
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Remark 2 When M = T*R™ is equipped with a canonical symplectic basis in 
cotangent coordinates z = (q, p) G T*R", the Hamiltonian vector field is given 
by Hamilton’s equations: 


■ dH ( f 

q = (q,p), 

<9p 


dH 

p = _ a5‘ <q,p) ‘ 


(4) 


Poincare discovered that the flow of any Hamiltonian vector field 5 is a 1- 
parameter subgroup of symplectic diffeomorphisms. Denoting such a flow by 
4>\j, this implies that for each fixed h £ R, (jffj is a symplectic map. 

Let Zj £ M be a point on the symplectic manifold and z (t) the integral 
curve to Xh such that zo = z(0). By definition of the flow, the mapping 

z(f + h)= 4> h H {z(f)) 


will propagate the solution from time t to time t + h. A symplectic algorithm 
with stepsize h, is the numerical approximation ifh of the Hamiltonian flow 
c bg : M —► M, which is an isometry of the symplectic form u> . Specifically, 
consider the exact solution z (t) of a Hamiltonian system for the time t £ [0, X 1 ], 
a discretization {U}^L 0 such that to = 0, ijv = T, h = T/N = t n + 1 — t n , and 
denote z n = z(t n ) for 0 < n < N. Let U C M be a convex open neighbourhood 
of z„ containing the target point z n+1 . 

With these hypotheses, we define a symplectic integrator as a map 


ifh-UcM ->■ U 

Z n 1 t z n+1 = ifh(z n ) 


smooth with respect to h and H , and such that ip’foJ = ui, where is the 
pullback of ifh defined by 

(^)z( v ,u ) =uj^ z) (Tif h (y),Tif h (u)), u, v G T Z M. (5) 

Remark 3 The vectors and Tip^fu) belong to the tangent space T^ h / Z ^M 

which, in general, is different from T Z M. Once we identify u z and u>^, h i z ) with 
wo and the tangent spaces T Z M and T^ h ( z \M with V, condition (5) becomes 

(v, Ju) = (Tip h (v),JTip h ( u)). 


or equivalently 

difh T T diph _ 
dz dz 

which is the well-known symplecticity condition for R 2 ™ viewed as symplectic 
vector space. It is worth nothing that in the last expression, z are already local 
coordinates. 

5 Poincare used the name of the “fundamental problem of dynamics”, which is to find the 
solutions of the “fundamental equations of dynamics in canonical form” [30]. 
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In an analogous way, we define an implicit symplectic integrator as a map 


ip h :U xU -»• U 

(Zra,Z n _|_l) I t Z n _(_i = (/P/j (z n , Z n _(_i) 

smooth with respect to h and H , and such that = oj. 

Since we are discretizing a flow, it is possible to consider an intermediate 
point z e U and two maps 6 7 tpi,ip2 : U —» U such that z = if>i(z n ) = ^(zn+i)- 
They let us rewrite the implicit scheme as iph{ z n , %n+i) = o ^i(z n ). 

The pullback becomes = (^>2 o = ip 1 0 (^2 )* 00 which produces 

the corresponding symplecticity condition in the tangent spaces by 



Note that on the left hand side of (6), J is an endomorphism on T Zn U and on 
the right hand side, J is an endomorphism on T Zrl+1 f7. 

Condition (6) says nothing about the mappings i/’i and ip2, but only that 
the composition is a symplectic map. The reader interested is referred to [16] 
for a deeper discussion on conditions imposed on ip\ and ip2- 

The discrete scheme iph is said of order r G N if, as h —> 0 

W<j>Hi Z n) - Vv»(zr i)ll = 0(h r+1 ), Z„ = Z (t n ) € U. 

There are several methods for constructing symplectic integrators that re¬ 
duce to finding symplectic maps between two different (closed) points on 
the integral curves of the vector field Xh ■ Here we are interested in the method 
of generating functions that we describe in the next section. 

3 Generating functions 

It seems that Poincare was the first author who systematically studied the pro¬ 
cess of obtaining symplectic maps by generating functions. However, the termi¬ 
nology was different in his work: Hamiltonian equations were called fundamental 
equations of dynamics in canonical form, symplectic maps were called canonical 
transformations'. Generating functions were the fundamental tool in his theory 
of integral invariants [30]. Poincare used these transformations to study periodic 
orbits of second genus in celestial mechanics. It may explain why this technique 
was not acknowledged immediatelly by the numerical community. From the 
numerical point of view, generating functions for symplectic integrators were 
systematically studied by Feng’s team in the mid ’80s [19, 20, 8, 9] and later, 
among many others, by Ge and co-workers in late ’80s and ’90s [13, 10, 11, 12]. 

6 In fact, they must be symplectic maps to have a consistent symplectic integrator [16]. 

7 Transformations which preserves the canonical form of the fundamental equations of dy¬ 
namics [30]. 
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Both the analytical and the numerical procedures use the same framework that 
we outline now. 

Let (Mi,wi) and (M 2 .^ 2 ) be two symplectic manifolds of the same dimen¬ 
sion. A map <j> : Mi —>• M 2 is called symplectic if (j>*u >2 = u>i, where, in general 
wi and u >2 are different symplectic forms. In our case, (Mi,wi) = (¥ 2 ,^ 2 ) are 
two copies of the same symplectic manifold, but we will preserve the subindices 
to keep record of what copy we refer to at every time. 

Consider the product manifold M = Mi x M 2 with canonical projections 
TTi : M —> Mi for * = 1,2, and define a two-form w e on M by 

UJq = 7T*Wi ^ 7T2W2 (7) 


We have the following results (see [1, sec 5.2] for the proofs): 

• (MjWq) is a symplectic manifold of dimension 4 n. 

• for any symplectic map (f> : M\ —» M 2 , the graph of <f>, denoted by and 
defined as 

= j(z,0(z)) G M I z e Mi,4>(z) e m 2 J , 

is a Lagrangian submanifold of M. This means We|r^ = 0 

• Since 9q is a closed 1-form by the identity d^elr^ = 0, using Poincare’s 
lemma, 9q is also an exact 1-form on Then, there exists a function 
S defined on the Lagrangian submanifold T^ such that its differential 
concides with the restriction of the 1-form 9q to P^ 

SM^—>-R, such that dS|r 0 = (8) 

S : T^, — > K is called a generating function for the symplectic map </>. 

Symplectic maps, generating functions and Lagrangian submanifolds are 
closely related. For instance, in a generic symplectic manifold M, any La¬ 
grangian submanifold A C M which is transverse to the fibers of the projection 
7 tm ■ T*M —> M can be locally parameterized by a suitable atlas of (local) func¬ 
tions Si for short times [14, 25, 4]. Since symplectic integrators are mappings 
close to the identity [h small), we are not concerned with the global behaviour 
of the Lagrangian submanifolds and all our analysis will be local. 

The standard procedure of the method of generating functions [19, 20] is as 
follows: 1) look for a suitable 1-form 9 G T*M such that d9 = w e ; 2) obtain the 
Lagrangian submanifold A C M, associated with 9 and a function S' : A —» K 
satisfying dS = 9 ; 3) solve the Hamilton-Jacobi equation on the Lagrangian 
submanifold; 4) design a numerical method for approximating such a solution 
giving the mapping z n+ i = iph(zn) for z n € M x and z n+1 € M 2 . 

Kang’s procedure for solving the Hamilton-Jacobi equation is to approximate 
the generating function using Taylor series expansions [19]. Menyuk used the 
Picard iteration to obtain such a function [26] and Channell and Scovel used 
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symbolic computing software [29]. As pointed out by Miesbach and Pesch in 
[27] , these approaches require higher order derivatives of the Hamiltonian which 
complicates the final scheme. 

To keep the notation simple, from now on lowercase variables belong to M\ 
and uppercase ones belong to M 2 , in particular z = z„ and Z = z n+ i. We 
use also Einstein notation in the definition of differential forms, i.e., ptdqi := 
£?=o Pi d Pi- 

3.1 The alternative method of Liouvillian forms 

Recently, the first author made some contributions on the subject [18], where he 
avoids the solution of the Hamilton-Jacobi equation and recovers a numerical 
algorithm from geometrical interpretation of the solutions. Indeed, solutions 
of the Hamilton-Jacobi equation on the product manifold M, belong to a La- 
grangian submanifold. This submanifold can be related to the characteristic 
bundle of the Hamiltonian vector field generated by a generalized generating 
function on M. The differential of this generalized function is actually a Liou¬ 
villian form, i.e. a 1-form A on M such that dX = u) q. A suitable projection 
7 r : M —> N to a 2n dimensional submanifold TV C M gives the right point 
where we must evaluate the discrete flow in order to have a symplectic inte¬ 
grator. Contrasting this construction with Kang’s procedure, we look for good 
values of the discrete flow on the Lagrangian submanifold defined by the Liou¬ 
villian form before the projection, instead of approximating the projection of 
the solution by Taylor series. Consequently, points 3) and 4) are not relevant 
in this construction. For more details the reader is referred to [18]. 

Our point of departure is the family of Liouvillian forms d a j 3,7 constructed 
in [18]. Consider two copies of the phase space Mi = T*Qi, M 2 = T*Q^ and 
define the product manifold M = Mi x M 2 equipped with the 2-form w@ given 
in (7). As we know, (M,Wq) is a symplectic manifold of dimension 4 n. Define 
local coordinates ( qt,Pi ) £ Mi and (Q i: Pi) € M 2 , * = 1, • • •, n, in each copy of 
the phase space and consider real numbers a, /3 ,7 £ E. With this notation, the 
family of primitive 1 -forms 0 a ,/ 3, 7 is given by 

d 11 . ;3. 7 — o:(pidqi -f- QidPi) (1 a)(PidQi qidpi) 

+/3(qidqi + QtdQi ) - q (pidpi + PidPi), (9) 

Since = <J# a ,/ 3, 7 for any values of the parameters, we can give a full set of 
3n different parameters {cti,/3j, 7 i }" =0 for every 1-form 0 q,,/ 3 )7 8 . Note that the 
elements associated with f3 and 7 belong to the kernel of the differential and they 
are known as the gauge elements of 0a 1 ( a, 7 . They do not modify the symplectic 
structure of (M, <J0 Qj( g )7 ) but the solution in the Lagrangian submanifold will 
be, in general, different for each combination of parameters. 

8 The case /?i = 7i = \/ai(l — on) with rv:,; £ [0,1] corresponds to the family 0^ constructed 
in [18] from the more elementary symplectic rotation on = cos 2 ( 0 ,), i = 1 ■ ■ ■, n, on M = 
Mi X M 2 . 
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Remark 4 Every point () £ R 3rt , * = 0, is associated to a 1- 

forrn whose differential is exactly the symplectic form . However, good values 
for a symplectic map approximating the Hamiltonian flow in a suitable way, 
belong to an open ball around the point = ( 1 / 2 , 0 , 0 ), associated to 

the mid-point symplectic map. In Section 4-1.2 we will find values for these 
parameters for a concrete example. 


All members of the family d a ,p n are 1-forms, locally closed on the graph 
r^cMofa generic symplectic map —>• (Aff,ujf). By Poincare’s 

lemma, there exists a generating function S = S a .g n , depending on (a, /?, 7 ) 
such that dS = 6 a: p^. The Lagrangian submanifold A aj( g i7 parameterized in 
local coordinates by the equation dS = 6 a ,f 3, 7 has explicit form 


dS_ 

dqi 

dS_ 

dpi 

dS 

dQi 

dS 

Off 


apt + f3qi 
-(1 -at)qi -'ypi 
-(1 -a)Pi+PQi 
aQi - 'yP l , 


( 10 ) 

( 11 ) 

( 12 ) 

(13) 


and must satisfy the homogeneous Hamilton-Jacobi equation H (A Qi ^ )7 ) = 0. 

Consider a 2n-dimensional submanifold of the product manifold N a _g^ C M 
with coordinates ( Qi , Pf), and dehne the projection 7 rg : M —> by 


n s = J o - n 2 ){A a ,g !j ). 


(14) 


The submanifold N a ,g rf is given in local coordinates by the equations 

Qi = aQi + (1 - a)qi + 7 (p t - Pi) 

Pi = api +(1-a)Pi +fi(qi-Qi). (15) 

Note that on the diagonal of M defined by 

A m — ^(.qi, Pi, Qi, Pi) ^ A1 | qi — Qi, Pi — Pi, Vf — 1 ■ • * , 

the submanifold N a _g rf is a standard symplectic submanifold with the induced 
symplectic form ua which fulfils Wq = n* s uj&. Moreover, given the inclusion 

i : M —»• M (16) 

(Qi,Pi) | -t ( Qi,Pi,qi,Pi ) (17) 


whose image is exactly the diagonal A^ = i(M, ui) we have ui = i*u> a- 
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3.2 The implicit symplectic integrators 

For (z, Z) £ M where z = ( qi,Pi ) £ Mi and Z = ( Qi,Pi ) £ M 2 , the family of 
symplectic integrators given in [18] is obtained by the following implicit scheme 9 

Z = z + h JVH(Tn s (z, Z)) (18) 


or in extended form 


Qi = qi + h— (aQi + a'qi + 7 (ft - P»), + a'Ps + /3(<7i - <?*)) (19) 

Opi 

dH 

Pi = Pi - h—(aQi + a'q i + 'y(p i - Pi),api + a'Pi +/3(qi - Qi)) 


where we used a 1 = 1 — a to simplify the expressions. We shall return to 
these expressions later, when we study the symmetric case. At this stage, it is 
worth to note that for (a, /?, 7 ) = (0,0, 0) and ( a , /3, 7 ) = (1, 0,0), we obtain the 
staggered Euler symplectic integrators A and B [7] 


dH 

Qi = Qi T (.Qii Pi) 

dpi 

dH 

Pi — Pi ^~r\ (.Qii P 0 
dqi 


dH 

Qi — Qi H"~ ( QiiPi ) 

dpi 

dH 

Pi — Pi ( Qii Pi) 

dqi 


( 20 ) 

( 21 ) 


The main remark on the choice of the method (18) is that projection tts given 
in (14) is a “rearrangement” of the data coded in the Lagrangian submanifold 
Ac*,/ 3, 7 in a close to symplectic submanifold N a ^^. This choice is rooted in the 
notion of geometrical solution to the Hamilton-Jacobi equation 

where the function H : M —> R of the Hamilton-Jacobi equation is evaluated 
on a Lagrangian submanifold A C M of a symplectic manifold (M, w) given by 
local equations 


A = < (q,p) £ M | q = q, p = 


dS_ 

<9q 


and differs from the energy Hamiltonian function, only by a constant (normally 
denoted by —E). For the non-evolutionary Hamilton-Jacobi equation, consid¬ 
ering the characteristic bundle, we evaluate the Hamiltonian vector field on the 
sub-bundle ( TX )" L C TM transversal to T A. An example is provided by the 
staggered Euler integrator of type A, with the difference approximation 


Z = z + hJVH(Q, p), (22) 

of the differential equation z = JVH( Q,p) associated to the Hamilton Jacobi 
equation H( q, P) = 0 for z = (q, p) and Z = (Q, P). 

9 Note that the projection (14) with local coordinates (15) is linear in (q t . p; ,Q r - Pj) and 
consequently it coincides with its linearization Tns- 
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4 Numerical examples 

All the members of the family (20) for which a € (0,1) produce implicit symplec- 
tic integrators that we can compute in an iterative predictor-corrector scheme. 

For the tests, we consider the perturbed Hamiltonian pendulum as toy ex¬ 
ample. The algorithm was implemented in python for testing the accuracy of 
the method and compared with the explicit method SABA 2 with coefficients 

(do = d 2 = 1/6, bo = bi = 1/2, ai = 2/3}, 

studied by Laskar and Robutel in [22] . We used the staggered symplectic Euler 
integrator as predictor and we compute, iteratively, the intermediate coordinates 
z = (Q,P) which are used to compute the value of z„ + i = (p„_|_i, q n +i) in k 
iterations. 

The general structure of the algorithm is the following 


Algorithm 1. 

1: 

! Setup the predictor 
zl°l =z n + hJSJH{ z„) 

2: 

for j = 0 : k do 

3: 

Icompute the point z = (Q, P) 

Q = aqM + (1 - a)q„ + 7(p n - p w ) 

4: 

P = ap„ + (1 - a)p w + /?(q n - q W ). 

5: 

Icompute the corrector 
ztl+d = z„ + hJVH{ z) 

6: 

end for 

7: 

Z„+i = z M 


4.1 The simple Hamiltonian pendulum 

The simple Hamiltonian pendulum obeys the Hamiltonian equation 

H{p,q) = ^p 2 -ecosq (23) 

where we have fixed e = 0.03. We performed several tests for different values of 
(a, /?, 7 ) and the step size h. The total number of steps N used in the simulations 
was variable but almost all the results used N = 10 5 . 

A full set of generic parameters (a, /?, 7) for our symplectic integrators in 
a Hamiltonian system with n degrees of freedom, has 3n elements. We are 
interested in some particular cases and the dimension of the space of param¬ 
eters is different for each one. Since it is not evident in a 1 degree of free¬ 
dom problem, we enumerate the dimension of the space of parameters for the 
cases used in the tests: 1 ) the case (a = 0.5,/? = 0,7 = 0 ) has null dimen¬ 
sion (is a point); 2) the case (a, /?, 7 = /?) has dimension 2 n; 3) the case 
{a,j3 = ■\Ja{\ — a), 7 = y/ot{l — a)) has dimension n; 4) the case (a = 0.5,/? = 
^/a(l — a), 7 = \Ja( 1 — a)) has null dimension, it is the point (0.5, 0.5,0.5); 5) 
the case (a = 0.5,/?, 7) has dimension 2 n. Below the reader will see that we 
look for optimal values of /? and 7 in an open ball centered at the origin of R 2rl . 
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Figure 1: A simple illustration of the natural error obtained by the symplectic 
mid-point rule. In this plot we assume known the vectors v = Xh(z) = z and 
V = X H (Z) = Z and we try to approximate V by computing X H ( 2 ^). Note 
that it coincides with Vi = ^{v+V) only if Xh is a linear vector field. Moreover, 
it coincides with V only when the solution z (t) is linear in t, or equivalently 
when Xh is constant. 


4.1.1 Preliminary tests 

We consider the classical symplectic mid-point rule as starting point since it 
is the algorithm corresponding to (a,/3, 7 ) = (0.5,0.,0.). The error associated 
to this case is easy to understand if we consider the projection of the orbit on 
the phase space (Figure 1). What we are looking for is some point z, as close 
as possible to the real orbit, such that on it, the real vector field is parallel 
to the numerical value. Of course, the point z coincides with the mid point 
i(z n + z„_|_i) when the Hamiltonian vector field is constant. 

When computing the symplectic mid-point rule with five iterations in the 
resolution of the implicit scheme, and comparing with the explicit integrator 
SABA 2 from [ 22 ], we found that the former is more accurate for the selected 
set of initial conditions (upper-right plot in Figure 4). However, it is significantly 
more costly because it is an iterative scheme. 

With this error reference’s framework, we consider the case where (a, /3, 7 ) 
are given by the simple symplectic rotation in [18]. The rotation corresponds 
to the parameters /3 = 7 = y/a(l — a) for a £ [0,1], The relevant observations 
follow: 

a) Variations in a from 0 to 1 result in a progressive deformation of the 
orbit which connects consistently the symplectic Euler methods A with B. The 
energy is well conserved with the well-known oscillations for Euler methods. The 
oscillations are of the same order of magnitude for every member of the family. 
In the symmetric case (a,/3, 7 ) = (0.5,0.5, 0.5), the oscillations are really large 
compared with the mid-point case (a,/?, 7 ) = (0.5,0.,0.) (Figure 2). 

b) Tests for stepsizes h between 0.001 and 0.1 exhibit a good behaviour and 
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Timestep h=0.1 le2 Timestep h=0.1 le2 


Figure 2: Oscillations around the energy integral of the numerical solutions 
computed with the new implicit method. Left column corresponds to initial 
condition (qo,po) = (0.8,0.0) and (qo,Po) = (3.1,0.0) for the right column. 
Different graphs in each panel correspond to the solution varying from a = 0.0 
to a = 0.5. Above: progresive deformation of the oscillations for the rotation 
0 = 7 = i/a(l — a). Below: the case 0 = 7 = 0.0. 


some deformation in the orbits is detected for 0.2 < h < 0.8 and /? = 7 = 
y/a(\ — a) (see right panel in Figure 3). For h > 1.0 the orbits experienced 
high oscillations and for initial conditions close to the hyperbolic fixed points 
the numerical solution went to the unbounded region. Deformation of the orbits 
in phase space corresponds to the oscillation of the numerical solution of the 
energy integral around the exact constant value of H. On the other hand, for 
small values of 0 and 7 the numerical solutions has a good behaviour (see left 
panel in Figure 3) and some tests for h = 10 are really satisfactory. 

The main remark from these preliminary tests is that parameters 0 and 7 
let us control the numerical error in the energy integral, and that values of the 
parameters close to (a,/3, 7 ) = (0.5, 0,0) have very small error. In the rest of 
the discussion we fixed a = 0.5 and we considered small values for the other 
parameters. This value for a is strategic since a controls the symmetry of 
oscillations around the constant energy for positive and negative directions in 
time. This is evident from the well-known fact that symplectic methods are 
reversible in time if and only if they are symmetric. Moreover, oscillations for 
a € ( 0 , 5 ) when 0 = 7 = 1 /a(l — a)) in positive time (a, t), are symmetric with 
respect to those for (1 — a, —t ) in negative time. 
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Figure 3: Orbits of the numerical solution for the perturbed Hamiltonian pen¬ 
dulum. We compare two different members of the family of implicit integra¬ 
tors (a,/3, 7 ) = (0.5,0,0) and (a,/ 3 , 7 ) = (0.1,0.3,0.3) vs. the exact solu¬ 
tion. There are four orbits in each plot, with initial conditions (qo,Po) G 
{(0.8,0.), (1.6,0.), (2.5,0.), (3.13,0.)}, stepsize h = 0.7 and k = 5. Note the 
deformation of the orbits at the right panel. 


4.1.2 The symmetric case a = 0.5 

We perform this analysis in two separate cases: the first one when 7 is free and 
f3 = 7 , and the second one when both /3 and 7 are independent. 

a) Let us fix a = 0.5 and /? = 7 , with 7 G [—0.5, 0.5]. In this case, the 
amplitude of the oscillations in the energy goes from a positive to a negative 
phase (upper-left panel in Figure 4). Moreover, for 7 ~ 0 and five iterations 
(for solving the implicit scheme) the amplitude of the oscillations is lower than 
oscillations of SABA 2 with a better behavior (lower-left panel in Figure 4). This 
fact is remarkable since our implicit symplectic integrator is a one step method. 
Moreover the tests on the behaviour of the error with respect to the stepsize 
show a better accuracy than SABA 2 (lower-right plot in Figure 4). 

A finer analysis for different orbits with initial conditions close to one of 
the hyperbolic fixed points gives us new insight for the understanding of our 
numerical scheme. 

b) Fix a = 0.5 and consider independent values for /3 ,7 G [—0.5,0.5]. We 
observe in this case, that each one of the parameters controls a different part of 
the oscillation around the energy integral. To show that, we choose an initial 
condition close to one hyperbolic fixed point. In an oscillation of a numerical 
solution given by a symmetric symplectic method (for instance SABA 2 ), three 
cusps arise 10 . In our results, /3 controls the central cusp in a very clear and 
definite way. This property is shown in the upper row of Figure 5. In this plot 

10 We need to check if other methods and other Hamiltonian problems behave in a similar 
way. 
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Figure 4: Comparison of the oscillations of the numerical solution around the 
constant energy integral, between the new integrator with a = 0.5 and the 
explicit SABA 2 integrator. All plots have initial condition (qo,Po) = (0.8,0.0), 
the stepsize in the three plot of energy integrals is h = 0 . 2 . Above-left: variation 
of the j3 = 7 parameters in {—0.5, —0.25, 0,0.25,0.5} Above-right: new method 
for /3 = 7 = 0 wompared to SABA 2 - Below-left: minimal oscillation around 
the energy integral of the new method for /? = 7 = —3.4 x 10 -5 . Below-right: 
generic members of the family have order 1. The dotted line corresponds to h 2 . 


we fixed 7 = —00.7 and we modified /3 visually such that the cusp be close to 
the constant energy. Then, we fix /3 and modify 7 , which controls the two cusps 
until we arrive at a very flat solution as we can see in the lower panel in Figure 

5. 

By a continuity argument we search for optimal values of /? and 7 , reducing 
the error in the numerical solution. This is done as a manual and visual process, 
computing the maximum variation of the error within an open square around 
(/3, 7 ) = (0, 0). We obtain a well defined region where such values belong. Max¬ 
imum variation in the energy integral for the perturbed Hamiltonian pendulum, 
with initial conditions (qo,Po) = (3.1,0.0), stepsize h = 0.1, perturbing term 
e = 0.03 and k = 5 iterations went to err < 5 x 10 -14 , which is remarkable for 
a one step integrator. Figure 6 shows the maximum variation of the error for a 
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Figure 5: Search of optimal values for f3 and 7 with a = 0.5. Above: variation of 
the P parameter changes the position of the central cusp. Below-left: variation 
of the 7 parameter changes the position of the sourounding cusps. Below-right: 
an optimal combination of parameters P and 7 reducing the error oscillations 
around err < 10 -13 . See left panel in Figure 6 for a zoomed version of this 
plot. The fixed oscillating orbit corresponds to the numerical solution of the 
SABA 2 integrator. All plots have initial condition (qo,po) = (3.1,0.0) and 
stepsize h = 0.1. 


subset of R 2 and the computed numerical energy compared with SABA 2 ■ The 
vertical line in the middle of the left panel in Figure 6 is the central cusp of the 
SABA 2 integrator. It is a zoomed version of the lower-right panel from Figure 

5. 

4.1.3 Interpretation of the parameters P and 7 

The fact that the new parameters P and 7 control the numerical error is an 
important result. However, the way we obtained these values is quite heuristic. 
In order for these parameters to be really advantageous, their behaviour with 
respect to the variations in the stepsize h must be investigated. 

To check their dependency on the stepsize h, we performed several tests 
with the same initial conditions as in the previous tests. We searched for the 
optimal values of P and 7 for h £ {0.05, 0.1, 0.5,1.0} and we obtained an almost 
perfect linear relationship p = h-b and 7 = h ■ c for b = —2.4978136594 x 10~ 3 
and c = —8.33321735568 x 10~ 2 (see upper row and lower-right panel in Figure 
7). To verify the linear dependence, we extrapolated and interpolated several 
values of the parameters for different h. The linear relationship holds in a 
very accurate way. The lower-right panel in Figure 7 is computed with the 
extrapolated parameters for h = 0.007. Figure 8 shows the values found by 
visual inspection for h £ {0.05,0.1, 0.5,1.0} and the line joining them. 
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Figure 6 : Search of the optimal values for (3 and 7 . Left panel shows the 
oscillations of the error in the energy integral computed with values for /?, 7 
found by the plot in the right panel. Vertical curves correspond to the numerical 
solution computed with SABA 2 - This is a zoomed version of the lower-right 
panel in Figure 5. Right panel shows the distribution of the maximal variation 
of the error in a small region of the /3 — 7 plane. The behaviour is equivalent for 
bigger and smaller regions: the minimal variation converges to a very particular 
region with radius around 10 ~ 13 . 


This fact reveals an important property giving a clue on the way a numerical 
symplectic integrator can preserve the energy (and any other) integral in a very 
accurate way. The extra elements in (Q,P) from (15) can be written as an 
approximation of the gradient vector field V z U(z) = — Jz considering 

7 (Pi ~ Pi) = 

P(Qi ~ Qi) = -hp (^) 


as the approximation hA{—J ^y 5 ) where 


A = 


l^n On 

0 „ -PIn ) ’ 


On, In G M nxn (M.). 


With this notation, the symplectic map (18) is redefined in terms of the 
numerical approximation of the gradient V H which belongs to the characteristic 
line bundle of the Hamiltonian flow 


Z 

h 


z 


JV Z H 


(T 


+ hA 



(24) 


In this form, the linear dependency of 7 and /3 on h has a geometrical mean¬ 
ing and provides some insight into the development of an analytical expression 
for the matrix A, which gives us the exact flow of the Hamiltonian system. This 
is the subject of a companion article [15]. 
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Figure 7: Energy integral of the numerical solutions of the new integrator 
compared to the explicit SABA 2 integrator with initial conditions (qo,po) = 
(3.1,0.0) and several stepsizes. For each value in the stepsize h, the values in 
the parameters /3 and 7 which minimizes the oscillations in the energy integral 
are proportional to h. Values of those parameters for the down-right panel were 
extrapolated by the constants of proportionality obtained from the other panels. 

5 Conclusions and perspectives 

In this paper we have implemented and tested the implicit symplectic integrators 
constructed in [18]. This implementation considers a simple and coarse iterative 
process which can be refined. The main goal here was to test the accuracy of 
the symplectic scheme. The construction process developed in [18] opens a 
new point of view on generating functions and symplectic integrators owing to 
several new results. 

i) First, we showed that gauge elements in the primitive 1-form 0 a ^enable 
the control of the numerical error in the energy integral and other integrals by 
construction. The set of parameters {a, /3, 7 } = {ai, /3j, 7 *}, i = 1, ■ ■ • ,n with 3 n 
elements, gives the set of all possible implicit symplectic integrators passing by 
z and Z, which are consistent with the projection (14) introduced in [18]. In the 
numerical tests, we searched for values of the parameters that produce the more 
accurate solution for the energy integral and the results are very promising. 

11 ) Second, numerical experiments support the existence of exact numerical 
sympectic integrators, contrary to Ge [11] and Ge and Marsden’s [13] claims. 
This subject is further elaborated in a short note [15] where Ge-Marsden’s claim 
is considered in the case of an implicit symplectic integrator. Their claim of the 
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Figure 8 : Left panel: The linear behaviour of the optimal values for f3 and 7 as 
a function of the stepsize h. Right panel, the optimal values of /3 and 7 for each 
h. increase the order of the integrator up to order four. 


non-existence of symplectic integrators exactly preserving the energy integral, 
relies on explicit symplectic maps generating the integrator. They also assume 
an extrapolation by Taylor series expansion. Our symplectic integrator is im¬ 
plicit, and the approximation is done by looking for internal points on the line’s 
flow. 

m ) Third, in the perturbed Hamiltonian pendulum, the parameters f3 and 7 
and the stepsize h satisfy a linear relationship. This fact, together with expres¬ 
sion (24) gives a clue on the way the difference equation recovers information 
from the Hamiltonian formalism to produce a very accurate numerical solu¬ 
tion. This subject is developed in [15] with a finer analysis on the symplectic 
transformation which better approximates the flow of a Hamiltonian vector field 
Xh■ Other Liouvillian forms as the one obtained from the Poincare’s generating 
function are studied in [17]. 

Further work is necessary to understand the subject of symplectic integrators 
using Liouvillian forms. However this series of papers gives a new point of view 
on the subject and on the control of numerical errors in symplectic integrators. 
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